Decreasing excitation gap in Andreev billiards 
by disorder scattering 
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We investigate the distribution of the lowest-lying energy states in a disordered Andreev billiard 
by solving the Bogoliubov-de Gennes equation numerically. Contrary to conventional predictions 
we find a decrease rather than an increase of the excitation gap relative to its clean ballistic limit. 
We relate this finding to the eigenvalue spectrum of the Wigner-Smith time delay matrix between 
successive Andreev reflections. We show that the longest rather than the mean time delay determines 
the size of the excitation gap. With increasing disorder strength the values of the longest delay times 
increase, thereby, in turn, reducing the excitation gap. 

PACS numbers: 74.45.+C, 73.21. La, 05.45.MT 
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A small metallic grain connected to a superconductor, 
commonly referred to as "Andreev billiard" (AB) fea- 
tures very intriguing electron dynamics that has been the 
focus of numerous studies, both theoretical and experi- 
mental 0, S H 0, S, 0, 1, [iJi^jn|JT2JTl_[ii, m [H, 
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view see [29]). The energy spectrum in such a grain is 
strongly influenced by the process of Andreev reflection 
of quasi particles at the contact between the supercon- 
ductor and the normal metal. At this " SN-interface" an 
incoming electron with energy e (counted from the Fermi 
energy Ep) is back-reflected as a "hole" with energy — e 
and nearly opposite momentum [l], [3^]. Such Andreev 
reflections result in the coupling between electron and 
hole excitations in the normal metal resembling those of 
electron-hole states in superconductors. In particular, 
the density of states (DOS) near the Fermi edge (Ep) is 
reduced and may exhibit an "excitation gap" (E±). De- 
tails of this reduced DOS are determined by the dynam- 
ics in the Andreev billiard which, in turn, depends on the 
boundary geometry, the position of the SN interface and 
on the potential surface in the AB. The distance E\ of 
the first excited state in the grain ("billiard") from the 
Fermi level (set equal to zero in the following) marks the 
size of the excitation gap in the energy spectrum. While 
being much smaller than the bulk gap A of the supercon- 
ductor, E\ <C A, Ei may considerably exceed the mean 
level spacing S, i.e. the average energy distance between 
adjacent eigenstates, thus signalling the appearance of a 
gap. 

Qualitative insights into the origin of spectral features 
in an Andreev structure, in particular the appearance of 
a gap, can be gained from a scmiclassical analysis. The 
semiclassical Bohr-Sommerfeld (BS) approximation [Tol 



111 . Il5l Il8l . |27| allows to relate Andreev-refiected peri- 
odic orbits with the energy levels of the Andreev bil- 
liard [see Fig. HJa)] . An energy eigenstate corresponds 
to a standing wave along such a periodic electron-hole 




FIG. 1: (Color online) (a) Geometry of an Andreev billiard 
(AB) consisting of a rectangular normal (N) conductor of di- 
mension (1.5W, 3W), where W is the width of the junction 
with the superconductor (S). The dotted (dashed) line depicts 
the electron (hole) part of a periodic electron-hole orbit cre- 
ated by Andreev reflection at the SN-interface. (b) A sample 
realization of the landscape of the disorder potential inside N. 



orbit with the action difference between electron and 
hole being quantized. This simple picture implies that 
the lowest energy E\ in the AB (i.e. the excitation 
gap) will be inversely proportional to the length of the 
longest electron-hole trajectories. However instructive 
the BS approach may be, it suffers from serious limita- 
tions, resulting from the assumption of strictly retracing 
electron- hole orbits. Corrections are due to short-range 
scattering off disorder 




quantum diffrac- 
tion [111, [16|, [191, [22J, [23J, [2JJ, deviations of the Andreev 
reflection from perfect backscattering 2^, 28| as well as 
due to contributions from trajectories that do not cou- 
ple to the SN-interface These corrections may 
turn out to be so large as to render a prediction for the 
excitation gap based on the BS approach unreliable. For 
example, the formation of a sizeable excitation gap in 
chaotic Andreev billiards as predicted by Random Ma- 
trix Theory (RMT) and verified numerically [nj, is not 
reproduced by the BS approximation [13J. 

In a realistic metal sample brought into contact with 
a superconductor, deviations from the ballistic limit by 
disorder scattering play an important role. If the elas- 
tic mean free scattering path £ is smaller than the lin- 
ear dimension of the metal grain, the trajectory between 
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FIG. 2: (Color online) (a) Disorder-averaged state count- 
ing function, (N(e))(, for four different disorder strengths 
V /Ep = 0.007, 0.03, 0.1 and 0.24 (colored solid lines) and 
Weyl estimate (black dashed line). The two lowest energy 
eigenvalues of the disorder-free system are marked by vertical 
bars, (b) Evolution of the mean gap {Ei)% (red triangles) and 
the root mean square deviation (blue squares) as a function of 
disorder strength Vo (in units of Ef). Horizontal lines mark 
the RMT predictions, (c) Statistical distribution of the low- 
est eigenvalue E\ for four disorder strengths Vo (colored lines) 
compared with the RMT distribution (black line) . (d) Depen- 
dence of the Wigner-Smith delay times on disorder strength 
Vo- Both the mean delay time, i.e., the dwell time (ra)? (blue 
squares), and the maximum delay times (r max )j (red trian- 
gles) are shown. The black dashed line shows the estimate 
(T d ) e = 2tt/NS' from H2]. 



two successive Andreev reflections at the SN-interface is 
dominated by disorder scattering in the interior of the 
grain rather than by ballistic scattering off the grain 
boundaries, ft has been suggested that the shortening of 
electron-hole orbits or, equivalently, of the average dwell 
time t<j between successive Andreev reflections by disor- 
der scattering would lead to an increase of the excitation 
gap as compared to that of a clean SN junction^, [l4j ]. 
Such a trend would qualitatively be in line with recent 
investigations [23j which have found that the gap in the 
ensemble averaged density of states of an AB increases as 
the mean free path decreases with respect to the clean, 
ballistic limit (under the assumption of constant average 
dwell times r<j). 

In this letter, we present numerical ab initio sim- 
ulations for a two-dimensional AB with disorder, em- 
ploying the Modular Recursive Green's function Method 
(MRGM) in combination with a wave-function matching 
technique at the SN interface 27|, l3l| . Disorder is rep- 
resented by elastic scattering off a potential distribution 
with short-range disorder with a correlation length ly 
small compared to the Fermi wavelength, ly/Xp = 0.12. 
Decohering processes such as inelastic scattering are ne- 
glected in the following. 



We choose a rectangular normal(N)-conducting cavity 
with dimensions (1.5W, 3W) where W is the width of 
the superconducting lead [see Fig. DJa)]. We construct 
the disorder potential [Fig. [IJb)] by decomposing the N 
region into two quadratic modules of dimension (1.5VY, 
1.5W) within each of which we employ a separable ran- 
dom potential, Vt(x,y) = Vt x (x) + V Sy (y) [£ x ,£ y denote 
two different statistical samples, jointly refereed to as 
£ = {£x)^y}]- This "trick" is employed for reasons of nu- 



merical efficiency, in particular for small Af 31| . We en- 
sure truly random scattering by destroying any unwanted 
separability by rotating by 180° the random potentials 
in the two squares relative to each other [see Fig. |TJ(b)] . 
The spatial correlation of the random potential is char- 
acterized by (Vt(x,y)V$(x + a,y)) x , y = (V^(x,y)V i (x,y + 
a))x,y — Vq x ex P( — a IW), with (• • • ) XjV indicating a spa- 
tial average over the whole disorder area and ly the cor- 
relation length. For a given realization £ the potential 
has zero spatial average, (V^(x,y)) x ,y = 0, and an ampli- 
tude, Vo = ([V^(x, y)] 2 ) x .y , which is chosen to be small 
compared to the Fermi energy, Vo/Ep < 0.2. For Vo —> 
the dynamics in the normal conducting part of the AB 
is entirely ballistic (no disorder scattering) and regular 
(due to the rectangular confinement). Calculations are 
performed with N = 24 open transverse modes fitting 
in the lead width W of the superconductor. The SN- 
interface itself is assumed to be ideal (no tunnel barrier) 
and the superconducting coherence length is assumed to 
be smaller than any other length parameters in the cav- 
ity. The superconducting gap A was chosen as 0.2-Ep to 
ensure that the energy E\ of the lowest-lying eigenstate 
fulfills E 1 < A. 

In the fully ballistic limit, i.e., in the absence of any 
disorder, Vq — 0, we find the lowest energy E\ to be 
four times larger than the AB's mean level spacing S [see 
Fig. [2a) for the lowest eigenenergies] . To investigate 
the influence of the disorder on the energy spectrum we 
now gradually increase the disorder amplitude Vo- For 
each value of Vo we calculate the full energy spectrum 
(below the superconducting gap A) for 500 different dis- 
order realizations £, and determine the ensemble aver- 
aged state-counting function N(e) (i.e., the integrated 
DOS) in the ensemble-average. For very weak disorder 
strength, V /E F = 0.007, we find (N(e)) 5 to be still very 
close to the fully ballistic limit of Vo — where the spec- 
tral density close to Ep is strongly suppressed relative to 
the Weyl estimate N(e) = ps for the DOS per unit area 
p = m c ff/(irh 2 ) [see Fig. HI a)]. Increasing the disorder 
amplitude Vo, however, gradually shifts N(e) towards the 
Weyl distribution [see Fig.[5{a,b)]. In particular, we find 
the size of the excitation gap to be reduced with increas- 
ing values of Vo, rather than increased. The reduction of 
Ei is a sizeable (factor 2 in the range < Vo/E-p < 0.2) 
and robust effect. For comparison we also show the gap 
as predicted by RMT for chaotic systems [see Fig. |2Jb)]. 
These RMT estimates are based on a numerical calcu- 



lation representing the internal dynamics of the normal 
conductor in the AB by an ensemble of 8000 symmet- 



ric random matrices of size M x M [lTf. M is assumed 
to scale with the ratio of cavity circumference C to the 
size of the SN junction W, M = N x C/W (for N = 24 
modes in the SN-interface we obtain M — 216). While 
the RMT value of the gap, E\ w 3.28(5, is in reason- 
able agreement with our numerical data for finite dis- 
order strength Vb ^ 0, significant discrepancies appear 
for the second moment (i.e. the variance) of the distri- 
bution y/ (-Ef)f [see Fig. H£b,c)]. The full quantum cal- 
culation shows first a steep increase in the variance with 
increasing disorder strength before levelling off, whereas 
the RMT result underestimates the width of the distribu- 
tion drastically. It should be noted that, strictly speak- 
ing, the limit of universality is only expected to hold for 
M > N > N 1 / 3 » 1 [10]. The latter limit is difficult 
to reach in any realistic simulation for a two-dimensional 
cavity. The fact that both the gap size and the variance 
stay at an almost constant value in a whole interval of the 
disorder strength, 0.1 < V$/Ep < 0.2, possibly points to 
a saturation effect due to the disorder-induced random- 
ization of otherwise boundary-specific scattering dynam- 
ics. The reduction of gap size and variance for stronger 
disorder, Vo/Ep > 0.2, may be related to a transition 
from weakly disordered scattering to diffusive or local- 
ized dynamics. 

The strong reduction of the gap size with increas- 
ing disorder points to a mechanism qualitatively differ- 
ent from the behaviour of the mean dwell time (r<i)f, 
which is only negligibly affected by increasing disorder 
[see Fig. H£d)]. To uncover the underlying physics we 
employ a rigorous approach that allows to relate the en- 
ergy spectrum of a quantum system to the dwell time dis- 
tribution that does invoke neither any semiclassical ap- 
proximation nor random matrix assumptions. Key to our 
approach is the relation between the Wigner-Smith (WS) 
time delay matrix Q and the scattering matrix S[33ll34j|. 



Q(s) = -ihS\e)dS(e)/ds 



(1) 



Equation (fTJ), well-known for unbound quantum systems, 
can be applied to the (bound) spectrum of an AB since 
an eigenstate of the AB occurs at an energy e for which 
the determinant 



det[l + S , (e)S' t (-£)] = 0, 



(2) 



where S{e) is the scattering matrix of the open, normal- 
conducting cavity with the superconductor replaced by a 
normal conducting waveguide of identical width W. Ex- 
panding S{e) around the Fermi energy (e = 0) for small 
e yields 
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FIG. 3: (Color online) Smoothed distribution of (Ei,Ei) 
pairs for three different strengths of the disorder potential, (a) 
Vo/E F = 0.007, (b) 0.03, (c) 0.24 (500 realizations of disorder 
were used). Black crosses in the contour plot mark the mean 
value of the distribution. Perfect (Ei,Ei) correlation would 
correspond to non- vanishing density only along the diagonal 
(drawn in the contour plot as guide to the eye), (d) Product 
of the mean gap size {Ei)£ and the mean of the maximum 
Wigner-Smith delay time (r max )5 as a function of disorder 
strength. The constant value hn/2 [predicted by Eq. |[5j] is 
shown for comparison. 



and, in turn, the approximate quantization condition for 
Andreev states [21] : 



1 



cxp ( 2 T £T r 



0. 



(4) 



The Wigner-Smith delay times r„ are the eigenvalues 
of Q. They correspond to "sticking" times inside the 
normal-conducting cavity between entering and leaving 
the cavity through the opening. Since in an AB the 
opening is replaced by an SN junction, t„ measures the 
dwcll-time between two successive Andreev reflections. 

The values of r„ (with n < N) provide a basis- 
independent measure for the sticking time of "quantum 
trajectories" without invoking model-specific assump- 
tions or semiclassical approximations. The only limita- 
tion of Eq. ((4]) is the error of order 0(e 2 ) due to the 
Taylor expansion and approximate resummation of the 
unitary operator S(e)S^(— e) [Eq. ©]. Equation (g]) re- 
lates the energy spectrum at small e to the largest delay- 
time eigenvalues. In particular, the size of the excitation 
gap E\ is determined by the maximal delay time value 



r max , such that 



Ei 



Hit 



2r n 



= Ei. 



(5) 



S{e)S\-e) = t + 2-eQ 
n 



exp 



(3) 



The disorder-averaged maximum delay time, (r max )^, is, 
indeed, monotonically increasing with increasing disorder 
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strength V [Fig. [2(d)]. In turn, Eq. ||5j) suggests that the 
disorder- averaged gap (Ei)£ will be reduced. 

To probe for the correlation between maximum delay 
time (r max ) and the gap size (£4) hinted at by Eq. (O, we 
have performed a statistical analysis of the distribution 
of (Ei,T max ) pairs for 500 disorder realizations (Fig. [3]), 
converted to (E\,E\) pairs using Eq. For perfect 

correlation we should expect the histogram to feature 
non-zero bins only along the diagonal {E\ = E\). Devi- 
ations from a perfect correlation, resulting in part from 
the Taylor expansion Eq. ([3]), provide a measure for the 
accuracy of the estimate E\ as compared to the exact gap 
size E\. For small disorder strength (Vq/E f — 0.007) the 
correlation between E\ and E\ is, indeed, near-perfect 
and non-zero bins occur only in a very limited range of 
values Ei, E x [see Fig. [3(a)]. With increasing disorder 
strength [V /E F = 0.03, see Fig. [3(b)] the maximum in 
the distribution shifts to smaller values of E\ and Ei 
while remaining correlated near the diagonal. Both ob- 
servations underscore that increased disorder decreases 
the gap which is, indeed, correlated with the maximum 
WS time-delay eigenvalue. For much stronger disorder 
[Vo/E F = 0.24, see Fig. [3(c)], the {Ex, Ex) correlation 
is diminished as off-diagonal bins become more signifi- 
cantly populated. While, on average, the connection be- 
tween the disorder-induced reduction of the gap and the 
increase of the maximal delay time (r max }j still holds, see 
Fig. [3d) , for individual strong disorder realizations this 
picture breaks down and small (large) gap sizes may well 
occur for systems with small (large) values of r max . 

The present simulations allow to directly inspect the 
effect of disorder scattering on the wavefunction densi- 
ties in the particle and hole sectors. The latter pro- 
vide a microscopic picture of the decay of correlations 
between gap and maximum delay time. In the ballistic 
limit Vq = 0, the electron and hole wavefunction tend to 
closely mirror each other [Fig. 0(a)] in agreement with re- 
tracing electron and hole orbits between two Andreev re- 




FIG. 4: (Color online) Electron and hole probability densities 
of the lowest Andreev eigenstate at (a) zero disorder poten- 
tial, (b) finite disorder strength Vo/Ef = 0.15 (one disorder 
realization) . 
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FIG. 5: (Color online) Distribution of Wigner- Smith delay 
times, P(t„) (colored lines) for an ensemble average over 
500 disorder realizations £ and different disorder strength 
(Vq/Ef = 0.007, 0.03, 0.1, 0.24). For each £ and V , the 
delay times were evaluated at 135 different energies in an in- 
terval _Ef±0.1A. Increasing the disorder strength Vo amplifies 
the long time tail of P(t„) (main plot, logarithmic scale), but 
concurrently produces much shorter delay times (see top right 
inset). The disorder- free delay times are indicated by vertical 
bars, with the lowest values starting at r « 2L/vf (the time 
of flight across the cavity length L — 3W and back) and the 
largest value at r as 131//-!; in- 



flections. With increasing disorder the similarity between 
wave components in the electron and hole sheet gradually 
disappears [see Fig. 2(b)]. This observation supports the 
picture that for strong disorder the wave function of the 
lowest AB eigenstate is largely determined by disorder 
scattering in the interior rather than by Andreev reflec- 
tions at the SN interface. Accordingly, the dwell time 
between two Andreev reflections looses significance. 

It is now instructive to inquire into the origin of the 
discrepancy to those models suggesting that the presence 
of disorder induces an increase rather than a decrease of 
the gap in comparison to its clean, ballistic limit 0, [3] • 
The key here is the disparate behaviour of the maxi- 
mum, (r max ){ = (max n= i ) jv r n){ , and average dwell time, 
(rrf)j = (J2n T n/N}^, the latter of which enters d, [3]- 
For the system under study here, disorder scattering is 
obviously able to "delay" for long-lived scattering states 
the interval between two successive Andreev reflections. 
The presence of disorder not only increases (r max )^, but 
also reduces the minimal delay time (r m i n ){ (see Fig. [5]). 
As a consequence, the distribution of delay times P(r n ) 
becomes "stretched", while leaving the mean value (t^ 
almost unchanged. The fact that the average dwell time 
(Trf)^ stays almost unaffected by the increasing disorder 
[see Fig. [5(d)] is in agreement with a general relation^ 
between the averaged trace of the matrix Q [see Eq. (JT])], 
(tiQ)^, and the mean spacing 5' of resonant levels in a 
(normal-conducting) scattering system, (trQ)^ = 2ir/5' 



5 



or, equivalently, (r^ = tt/NS (note that 8' = 2 x S, 
with 5 being the mean level spacing in the correspond- 
ing AB). Consequently, the mean dwell time should be 
entirely independent of the disorder configuration. This 
is, indeed, very accurately confirmed by our numerical 
results for (r^ [see Fig. H£d)]. In turn, (t<j)^ is un- 
suitable for characterizing correlations between gap size 
and disorder strength, because of its independence of Vq. 
Therefore, relating the gap size to the mean dwell time 
(rd}{ also fails to account for the gap reduction observed 
here (35|. Clearly, the present results do not preclude 
an increase of the excitation gap with increasing disorder 
for particular boundary shapes, e.g. for a gapless exci- 
tation spectrum in the absence of disorder. The present 
analysis suggests, however, that also in such systems the 
behaviour of the longest WS time delay eigenvalue will 
control the behaviour of the gap. 

The results in Fig. \5\ demonstrate that for a disordered 
cavity the strength of the disorder (Vq) does have a cru- 
cial influence on the distribution of delay times P(t u ) 
(in particular for long times). For chaotic cavities it 
was found both classically [361 ] and quantum mechani- 
cally [l^, 31, 37 1 that the long time tail of delay times 
does not depend on certain system-specific parameters 
as, e.g., the Lyapunov exponent. We therefore expect 
that the statistical distribution of the gap size undergoes 
a crossover between the present case of a disordered AB 
and the case of a chaotic AB. It would be interesting to 
study such a crossover numerically, e.g., by tuning the 
disorder correlation length ly from the diffractive limit 
of ly <C Af to the ballistic (chaotic) limit of ly 3> Xf- 

With the help of Fig. [3 we can furthermore explain the 
loss of correlations among pairs (E\,Ei) for strong disor- 
der [Fig. [3(c)]: The amplification of the maximal proper 
delay times (r max )j by the increasing disorder translates 
into an increase of the high frequency components in the 
elements of the scattering matrix S(e). As, however, the 
estimate of the gap size, E\, relies in part on a Taylor 
expansion of 5(e) [see Eq. ([3])] which can only capture 
weakly energy dependent (i.e., low- frequency) compo- 
nents, the accuracy of E\ deteriorates with increasing dis- 
order strength, thereby gradually diminishing the corre- 
lations among pairs {E\ 1 E\). The behaviour of the mean 
values (E±)c and (E\)^ can be understood by considering 
the distribution of values ds = [E\—E\)/2 (correspond- 
ing to a projection of the distributions of Fig. (3(a-c) on 
an axis perpendicular to the diagonal E\ = E\). As we 
have verified numerically (not shown), the width of this 
distribution, vax(dE) = y (d%)^, increases with increas- 
ing Vq, while its mean value stays almost unaffected by 
the disorder strength at <C S. We speculate that the 
errors due to the Taylor expansion and the resummation 
of S^(e)S(— e) [see Eq. ([3])] are randomly distributed and 
thus cancel out in an average over many disorder realiza- 
tions. This would explain why the averaged values (Ei)^ 



and (Ei)^ = (T max )£ remain correlated [see Fig. [3(d)] 
while the correlation between individual pairs {E\,Ei) 
breaks down. 

In summary, we have numerically calculated the energy 
spectrum of electron-hole states in a rectangular Andreev 
billiard with a tunable disorder potential. In apparent 
contrast to qualitative models based on the mean cav- 
ity dwell time (rd)^, we find a decrease of the gap size 
when increasing the disorder amplitude. We show that 
this decrease is controlled by the disorder dependence of 
the largest Wigner-Smith delay time r max between subse- 
quent Andreev reflections at the SN-interface. The aver- 
age dwell time (Td)j, on the other hand, only depends on 
the mean level spacing, and is thus neither correlated to 
the evolution of the gap size nor to the disorder scatter- 
ing strength. Stronger disorder, however, drastically in- 
creases the value of T max for long-lived scattering states. 
For sufficiently strong disorder the correlation between 
the gap size and r max eventually breaks down for indi- 
vidual disorder realizations, as the eigenenergies of the 
system are then more strongly influenced by the specific 
disorder potential rather than by the Andreev reflection 
process. 
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